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ABSTRACT 

We present a radiative magnetohydro dynamics simulation of the formation of an 
Active Region on the solar surface. The simulation models the rise of a buoyant mag- 
netic flux bundle from a depth of 7.5 Mm in the convection zone up into the solar 
photosphere. The rise of the magnetic plasma in the convection zone is accompanied 
by predominantly horizontal expansion. Such an expansion leads to a scaling relation 
between the plasma density and the magnetic field strength such that B oc q 1 ! 2 . The 
emergence of magnetic flux into the photosphere appears as a complex magnetic pat- 
tern, which results from the interaction of the rising magnetic field with the turbulent 
convective flows. Small-scale magnetic elements at the surface first appear, followed 
by their gradual coalescence into larger magnetic concentrations, which eventually re- 
sults in the formation of a pair of opposite polarity spots. Although the mean flow 
pattern in the vicinity of the developing spots is directed radially outward, correlations 
between the magnetic field and velocity field fluctuations allow the spots to accumulate 
flux. Such correlations result from the Lorentz-force driven, counter-streaming motion 
of opposite-polarity fragments. The formation of the simulated Active Region is ac- 
companied by transient light bridges between umbrae and umbral dots. Together with 
recent sunspot modeling, this work highlights the common magnetoconvective origin of 
umbral dots, light bridges and penumbral filaments. 



1. Introduction 



Active regions (ARs) are the most prominent manifestation of the large scale solar magnetic 
field in the photosphere of the Sun. Understanding the underlying flux emergence process is a crucial 
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step toward understanding the connection between the dynamo processes in the solar convection 
zone and magnetic activity in the photosphere and above. 

Observations of active regions on the solar surface indi cate that, in term s of magnetic flux 
content, there is a continuous spectrum of active region sizes (lZwaanlll978l . 119871 : lHarvey &: Martin 



19731 : ISchrijver Zwaanl 120001 : lHagenaarl l200ll : lHagenaar et al 



200. 



At the small end of the 
spectrum, ephemeral ARs (which have a flux of 3 x 10 18 — 1 x 10 20 Mx) manifest themselves in the 
solar photosphere in term s of mixed polarity field and magnetic bright points but not necessarily 
pores (jCheung et al.ll2008l ). Small ARs, which have fluxes of 1 x 10 20 to 5 x 10 21 Mx in each polarity, 
contain solar pores. Large ARs, which have even more magnetic flux (up to 4 x 10 22 Mx), contain 
sunspots with penumbrae. 



Since the original suggestion by iParkerl (|1955l ). there has been a large body of theoretical 
and observational work supporting the scenario that sunspots (and ARs) form as a result of the 
buoyant rise of magnetic flux bu ndles from the solar con yection zone to the sola r atmosphere. On 
the observational side, the work of I Zwaanl (|l978l . ll987l ) and lStrous k, Zwaanl (|1999l ) showed that ARs 
form as a consequence of a succession of small-scale flux emergence events on the solar surface. The 
observational studies in these papers noted that although magnetic flux emerges as small bundles, 
these bundles subsequently migrate and coalesce in a systematic fashion, eventually leading to 
formation of pores and sunspots (i.e. p atches of predominan t ly one polarity). Recently, an analysis 
of spectropolarimetric observations by ISchlichenmaier et al. shows how the accumulative of 

flux by a solar pore resulted in the formation of penumbral filaments. 

Numerical magnetohydrodynamics (MHD) simulations of the rise of buoyant magnetic flux 
tubes have yielded a wealth of information regarding the flux emergence process. For example, 
numerical models of the rise of toroidal flux tubes from the bottom of the convection zone based 

are able to reproduce 



on the thin flux tube approximation (R oberts fe Webb 



1978; 



Spruit 



1981 



global pr operties of ARs such as th eir tilt angles ( D'Silva fe Choudhuril 1993 ). their emergence 



latitudes (Caligari et al 



ties (IFan et al 



1993 



1995, 



1998 ) and asymmetries between the leading and following polar! 



Moreno-Insertis et al 



1994 



Caligari et al.l Il995l . Il998l ) . Due to the invalid 



ity of the thin flux tube approximation in the uppermost ten or so Mm of the solar convection 
zone, such studies have limited applicability for examining how ARs for m on the solar surface - 
Complementing t his line of work, mul t idimensional MHD simulations (e.g 



Isobe et al 



2006 



200 



5 



Shibata et al. 



1989; 



Archontis et al. 



2004; 



Forbes Sz Priest 



Manchester IV et al 



1984 



20041 ; iMagara 



Fang et al.ll2010l ) into the solar atmosphere have greatly improved our understanding of how 



the buoyant rise of magnetic fields into the solar atmosphere affect local dyna mics. For a his- 
torica l perspective of such wor k, th e reader is r eferred to the review articles by iMoreno-Insertis 
(|l997h . lFisher et al.1 (|200ol ;) andH (booi hood ). 



Recent models that include the treatment of convective motions and radiative heating/cooling 
(either by solving the radiative transfer equation or by parameterized terms in the e nergy equation) 
have shown these two effects to be very important to the properties of flux emergence (jStein Sz Nordlund 
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2006 :ICheung et alJl2 007. 20 081: llsobe et al.ll2008l: lAbbettll2007l: iMartinez-Svkora et alJIgOOa . 12003: 



Tortosa-Andreu fc Moreno-Insertis 



2009 



Yelles Chaouche et al.1 12009 : iFang et all l2010h . Despite 



these advances, a comprehensive numerical model of the formation of an AR resulting from flux 
emergence has been lacking. This is a consequence of the vast range of length and time scales 
encountered in the convection zone as consequence of a density variation by about six orders of 
magnitude from the base to the top. To cope with this difficulty numerical models focus either on 
the deep convection zone leaving out the upper most 10 - 20 Mm or on the upper most 10 Mm 
including the photosphere. 



While the former typ ically utilize the anelastic approximation (jAbbett et al.l 12000 : iFan et al 



iYP 

20031 ; iJouve Sz Brunll2009l ) the latter have to be based on fully compressible MHD. Regardless of the 
differences a common challenge in both cases is to explain the formation of rather coherent sunspots 
with several kG field strength from magnetic field which has risen through a highly stratified 
atmosphere and consequently should have weakened considerably through expansion. Addressing 
this aspect through a self-consistent numerical simulation capturing the rise of magnetic flux from 
the base of the convection zone into the photosphere is currently out of reach, however much 
can be learned from a numerical simulation containing the last 7.5 Mm beneath the photosphere. 
While the geometric extent is only 3.5% of the convection zone depth, the density drop of 3 orders of 
magnitude (corresponding to ten pressure scale-heights) is comparable to the drop in the remaining 
96.5% of the convection zone. In this paper we investigate the flux emergence process in the upper 
most 7.5 Mm of the convection zone resulting from the advection of a semitorus-shaped flux tube 
across the bottom boundary of our computational domain. 



2. Simulation Setup 



The radiative M HD simulation was carried out with the MURaM code (jVogler et al.ll2005l : 
Rempel et al.ll2009bl ). which takes into account radiative energy transport in the energy equation 
(assuming local thermodynamic equilibrium) and a realistic equation of state. The MURaM code 
has been used to model a variety of magnetic structures in the solar photosphere and convection 
zone. The Cartesian domain spans 92 Mm x 49 Mm in the two horizontal directions and 8.2 Mm in 
the vertical direction (with horizontal and vertical grid-spacing of 48 and 32 km respectively). The 
base of the photosphere (z = 0, the mean geometrical height where the Rosseland optical depth 
tr has a value of unity) is located 7.5 Mm above the bottom boundary of the simulation domain. 
Before the introduction of a magnetic flux tube, a purely hydrodynamic simulation was performed 
to allow the radiatively driven convection to relax to statistical equilibrium. 

To mimic the rise of magnetic flux into the topmost layers of the convection zone, an axisym- 
metric, twisted flux tube with the shape of a semitorus is kinematicaily advected into the domain 
through the bottom boundary. Similar time-dependent boundary conditions to kinematicaily ad- 
vect magnetic field into the top layers of t he co nvecti on zone have prev iously been used for flux 
emergence simulations. IStein h Nordlundl ([2006) and IStein et al.1 pOlCj) used a boundary condi- 
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tion whi ch advect purely horizo n tal fie ld (of uniform field strength and orientation) through upflow 



regions. iMartmez-Svkora et al.l (|2008l ) used a boundary condition which advected an horizontal 
twisted flux tube. The boundary condition for the current simulation is implemented by specifying 
the time-dependent fluxes of vertical mass, momentum, energy and magnetic field (i.e. terms acted 
upon by the divergence operato r in the conservation form of the MHD equations) consistent with 
the rise of the torus. Following iFan Gibson! ([2003), the magnetic field distribution of the flux 
tube in a spherical coordinate system with the torus centered at the origin is given by 



B 



A 



V x [ 



A(r,6) 



r sin! 



+ BJr, 



— XaB t exp (-co 2 /a 2 ), and 
(r sin 6)~ 1 aBt exp (— uj 2 /a 2 ). 



(1) 

(2) 
(3) 



uj = (r 2 + R 2 — 2rRs\n6) 1 / 2 is the distance of a point r from the toroidal axis of the tube, and 
a and R are the semi-minor and -major radii of the to rus, respectively. Th e dimensionless twist 
parameter A is equivalent to the parameter a~ 1 q used bv lFan &: Gibson! ( 20031 ). For this simulation, 
we took A = 0.5, R = 16 Mm and a = 3.6 Mm. Bt = 94 kG such that the toroidal field strength at 
the tube axis is B^(uj = 0) = aR~ l B t = 21 kG (corresponding to a plasma j3 of 140). For uj > v2a, 
B tu b e is set to zero. The total toroidal flux content of the tube is <&o = 7.6 x 10 21 Mx. Averaged 
over the cross-section of the tube (i.e. < uj < V2a), the mean toroidal field strength is 9 kG. This 
is somewhat stronger, but still comparable to field strengths of a few kG expected from thin flux 



tube simulations that begin with 100 kG at the base of the convection zone (see, e.g.. lCaligari et al 
19951 ). Thin flux tube simulations also predict a rise speed of between 0.5 and 1 km s _1 at this 
depth. For this simulation, we chose to impose a rise speed of 1 km s" 1 . 

The axis of symmetry of the torus is parallel to the y-axis of the Cartesian domain so that the 
resulting magnetic loop in the domain has its axis lying in the x-z plane. At the bottom boundary, 
the gas pressure within the torus is given by p gas = (pbot) — B 2 /8ir, where (pbot) = 2.44 x 10 9 dyne 
cm -2 is the mean gas pressure at the bottom boundary. The inflow specific entropy s of the plasma 
in the torus is set to the same value as that of ambient convective upflows. This choice of the 



thermodynamic properties results in a relative density deficit of A^/(^>bot) 
(71 = §1^1 a is the first adiabatic exponent) at the tube axis. 



(p-nY 



6 x 10" 



Beginning at t = 0, the magnetic torus is kinematically advected into the domain with a 
constant rise speed of 1 km s" 1 until the top half of the torus has traversed through the bottom 
boundary. Beyond this point in time (i.e., t > 5.9 hrs), the velocity at the bottom boundary within 
the tube is set to z ero. Outside o f the tube, the bottom boundary condition allows for smooth inflows 
and outflows (see IVogler et al.l 2005 for details). Mass flows across the top boundary are allowed 
(upflows are unimpeded and downflows are damped). By virtue of the low densities near the top 
boundary, the mass flux across the top has a negligible effect on the mass content in the simulation 
domain. The magnetic field at the top boundary is matched to a potential field. Periodic boundary 
conditions apply at the side boundaries. Following iRempel. Schiissler. &; Knolkerl ()2009bl ). we limit 
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the strength of the Lorentz force in low-/? regions (/3 < 0.05) to limit the Alfven speed to a maximum 
of 31 km s _1 . This artificial limiting of the Lorentz force mainly acts within the central umbral 
regions of sunspots in mid-to-high photosphere, where local Alfven speeds reach up to 10 3 km s _1 
and would impose prohibitively small time-steps on the explicit scheme of the code. 



Recently. iMacTaggart &: Hoodl (12009) modeled the rise of buoyant toroidal magnetic flux tubes. 
Their simulation setup is different from ours in a number of ways. First of all, they used a plane- 
parallel background convection zone and atmospheric model to mimic the stratification in the Sun 
but do not include convective flows nor radiative transfer. Secondly, the toroidal flux tubes in 
their simulations were embedded as buoyant, stationary structures as an initial condition. An 
interesting result from their calculations is that the choice of a toroidal flux tube (as opposed 
to an ori ginally horizontal tub e ) faci litated the emergence of the axis of the tube into the pho- 
tosphere. IMacTaggart &: Hoodl (|2009l ) attributed this to the ability (in the toroidal tube case) of 
plasma to drain down the legs of of the tube to enhance buoyant at the tube apex. In section 13.31 
we showed that convective flows acting on the rising magnetic plasma have a similar effect. 



Simulation Results 



3.1. Photospheric Emergence and Active Region Formation 



Before proceeding to discuss the physics underlying the flux emergence and AR formation 
process, we describe how the emerging flux region evolves in time. Magnetic flux begins to emerge 
into the photosphere beginning at t = 2 hrs and emergence progresses for several hours. A time 
sequence of the grey intensity and maps of the vertical field (sampled at tr = 0.1) are shown in 
Fig. [TJ In the intensity map at t = 3 hr, elongated granules begin to appear. This stretching 
of the convective cells is due to the horizontal expansion of the rising magnetic plasma. A more 
detailed discussion of the physics of horizontal expansion is given in section 13.21 Due to the 
subsurface horizontal expansion and undulation of the field lines by convective flows, the emerging 
flux region covers an extend surface area within which small-scale bipoles emerge with a systematic 
orientation (see also ICheung et a l. 2008) . By t = 6.7 hr, small pores begin to appear in multiple 
locations in the intensity map (c.f. IStein et al.l 120101 ) . A corresponding synthetic 'magnetogram' 
(B z sampled at r = 0.1) is shown in Fig. [2j Two hours later, two spots have appeared near 
the horizontal positions coincident with where the subsurface footpoints of the initial semitorus 
are rooted. As time progresses, the intensity maps show the spots growing in size. A synthetic 
magnetogram at a later time shows two large, vertical concentrations of opposite polarity field 
at the locations of the spot (see Fig. |3|). Throughout the simul ation, the two spot s never attain 
fully developed penumbrae which is consistent with observations (jZwaanlll978l . Il987l ) since the two 
large magnetic concentrations have (over the course of the simulation) at most 4 — 5 x 10 21 Mx. 
Nevertheless, transient penumbral filaments, umbral dots and even light bridges are found in the 
synthetic intensity maps. 
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3.2. Subsurface Rise and Expansion of the Magnetic Tube 



Figure H] shows the structure of the torus as it begins to erupt from the convection zone into 
the overlying photosphere. The subsurface roots of the torus at the bottom of the computational 
domain are shown in the magnetic map at a depth of 7.5 Mm (left panel). Field lines traced from 
the opposite polarity roots diverge upwards to covered an extended horizontal region. The right 
panel shows the associated emerging flux region at the surface, where the magnetic field lines have 



(see e.g. 


Pariat et al. 


2004; 


Cheune et al.l 


2007; 


Centeno et al. 


2007; 


Martinez Gonzalez et al. 


Martinez Gonzalez & Bellot Rubio 


2009; 


Cheung et al. 


2008; 


Ishikawa et al. 


2oioh. 



The number of pressure scale-heights spanned between z = —7.5 Mm and z = (i.e. the 
photospheric base) is N(H P ) = 10. To reach pressure equilibrium with the surroundings, the rising 
magnetic structure expands strongly in the horizontal directions. This is an important stage in 
the evolution for reasons of hydrostatic balance and mass conservation. Since the rise time of the 
flux tube is much longer than both the free-fall and sound-crossing times (both ~ 10 2 s) across 
the layers of the convection zone captured in this simulation, its rise does not significantly alter 
the mean stratification. On the other hand, the injection of mass from the rise of the half-torus 
is equal to 40% of the original mass in the entire domain, which is roughly equal to the total 
mass content of the top 6 Mm of the convection zone (as captured in the computational domain). 
Thus it is not surprising that the rising magnetic plasma displaces a substantial fraction of the 
original mass in the domain and fills the near-surface layers during the first stage of the active 
region formation process. As a consequence, the top of the magnetic torus takes on a flattened, 
pancake-like structure in near- surface layers (see Fig . [5l). A similar behavior is also reported in the 
flux emergence simulations by lToriumi Yokoyamal (|2010l ). 



Given the field strength of the original magnetic torus at a depth of 7.5 Mm, what strengths 
can we expect for the field emerging into the photosphere? This question can be addressed by 
considering a number of different simplified scenarios. For a purely horizontal magnetic flux tube 
expanding in directions transverse to its axis, conservation of mass and conservation of magnetic 
flux give the linear scaling relation B/g = constant. A different scaling relation can be derived by 
considering how an upwelling in the stratified convection zone modifies horizontal field threading 
the rising plasma. Let the initial horizontal field be B = Bx and for simplicity assume that 
the upflow is axisymmetric about the (vertical) z-axis and centered at the origin (the absolute z 
coordinate is unimportant). Let the rates of expansion in the horizontal and vertical directions be 
a = dv x /dx = dv y /dy and ea = dv z /dz respectively. Here e is a measure of the flow anisotropy 
(e = 1 for isotropic expansion). From the ideal Induction Equation, the rate of change of the 
magnetic field attached to a Lagrangian fluid element is 



DB 

~Dt 



-(V-v)B + (B- V)v, 



(4) 
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which in this case reduces to 

___ = _(l + e ) a . (5) 

The corresponding Lagrangian form of the continuity equation is 

Ding 



Dt 



V • v (6) 



= -(2 + e)a. (7) 
Combining equations © and ([7]), we obtain the scaling relation 

B(X g(l+e)/{2+e). (8) 

For isotropic expansion (e = 1), -Boc g 2 /' 3 . In the case where the horizontal expansion rate is much 
larger than that of the vertical expansion rate we have £ < 1, so that B oc ^J~q. 

For a field strength of 10 kG at g = 4 x 10 -4 g cm -3 (density at z = —7.5 Mm), the linear 
scaling relation gives, for g p ^ ~ 4 x 10 g cm -3 , a photospheric field strength of ~ 10 G. This is 
evidently too weak compared to what is found in the simulation and to what is observed, both of 
which give emerg ing horizontal field strengths on the order of a few hundred G (jLites et al.lll998l : 



Kubo et al J 12003 ). 



l + E 

The second scaling relation B oc g 2 +? yields more consistent values for emerging horizontal 
field strengths. For e = (zero expansion in the vertical direction) and e = 1 (isotropic expansion), 
a density drop of 10 3 yields horizontal photospheric field strengths of B = 100 G and B = 300 
G, respectively. Figure [6] shows a joint probability density function (JPDF) of the horizontal field 
strength versus mass density sampled at the mid-plane x = between t = 1 and t = 8 hours. Points 
with low values of specific entropy (s < 9.75 x 10 s ) have been excluded since they correspond to 
plasma that has undergone radiative cooling at the surface. The maximum magnetic field strength 
of 21 kG corresponds to the field strength at the axis of the torus that is introduced through the 
bottom boundary. The solid and dashed lines correspond to the scaling relations B oc ^fg and 
B oc £> 2 / 3 respectively. The latter scaling relation is clearly too steep and underestimates the field 
strengths at photospheric densities (g p h k 4 x 10~ 7 g cm -3 ). The square-root scaling relation 
provides a much better match to the JPDF. This result is consistent with the fact that the plasma 
experiences predominantly horizontal expansion (e close to zero) during its rise to the surface (see 
Fig. EJ. 

Figure [7] shows how the arrival of buoyant, magnetic plasma at the photosphere leads to a 
transient pressure excess which drives diverging horizontal flows about the flux emergence site. 
The coincidence (spatially and temporally) of this pressure enhancement with the onset of large- 
scale diverging flows (beginning at about t = 1 — 2 hrs) is evidence that the flows are driven by the 
associated horizontal pressure gradient. AP/(P) is of the order 0.1 — 0.5, so that the horizontal 
flows have Mach numbers of similar amplitude. After the initial emergence phase (t > 6 hr), 
the divergent horizontal flows at the periphery of the magnetic complex persist, albeit with lower 
amplitude (\v x \ ~ 1 — 3 km s _1 ). This raises the question of how magnetic flux is still capable of 
accumulating in the two spots, which will be addressed in section l3~4l 
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3.3. Discharge of Mass from the Rising Magnetic Structure 



Given that the near-surface field is rather dispersed, the next stage of active region evolution 
is the coalescence of the dispersed field into coherent concentrations. In order for the fragments 
to collect together, excess mass must be unloaded from the emerged field lines. To illustrate the 
amount of mass that is eventually removed, we compare the mass carried within the original half- 
torus and the mass remaining within the two opposite polarity concentrations at t = 8 hrs. Let 
B(r, z) and ~g(r,z) be azimuthal averages of one such flux concentration about its vertical axis. 
Using these, we calculated surfaces of constant magnetic flux and the corresponding mass content 
contained within the flux surfaces. The mass contained within a volume defined by the flux surface 
<&(r, z) = 4 x 10 21 Mx (summed over both spots) at t = 8 hrs is only 12% of the mass contained 
in a sub- volume of the original half-torus with the same flux content. Thus most of the original 
injected mass has somehow escaped to the surroundings. 

A number of physical mechanisms may be invoked to explain the discharge of mass from the 
original magnetic structure. First of all, magnetic diffusion of the torus allows mass transfer across 
magnetic field lines. However, this would result in an increase of mass in the semitorus rather 
than a decrease and thus can be ruled out as the underlying mechanism for mass discharge. As a 
second possibility, we consider the mass flux through the top and bottom boundaries. By virtue 
of the low densities at the top boundary, the vertical mass flux there is negligible for the mass 
budget of the semitorus. At the bottom boundary, vertical (and horizontal) velocities within the 
torus {uj < \^2a) are set to zero after the initial prescribed rise. Downflows are permitted at other 
magnetic footpoints at the bottom boundary but since the majority of the magnetic flux at that 
layer remains within the torus, mass flux through the bottom boundary can also be ruled out as 
the main discharge mechanism. A third possibility is the removal of mass via outflows associated 
with the horizontal expansion (see section I3.2[) . While it is true that outflows carry mass away 
from the center of the emerging region, these flows also advect magnetic field with them and cause 
a weakening of the field. So outflows alone are insufficient to explain how mass is removed from 
the field lines. 

The responsible mass removal mechanism is illustrated in Fig. [HJ The left panel of this figure 
shows a schematic representation of convective flows undulati ng a field line which has one end 
anchored deep below the photosphere. As already reported in ICheung et al.l ()2008l ). expulsion of 
flux from the granular centers lead to the encounter and cancellation of opposite polarity field. Such 



oppo site polarity pairs are connected below the surface in the form of U-loops (see also IStein et al 



2010I ). The first consequence of this magnetic reconnection between the opposite polarities is 
the creation of an O -loop, which d ischarges mass from the original field line. This is akin to the 
process suggested by iParkerl (| 19941 ) for the discharge of mass from and intensification of magnetic 
fields at the bottom of the conve ction zone. Base d on observations of emerging flu x regions with 



Spec tro-Polarimeter instrument (ILites et al.l 120011) of t he So l ar Op tical Telescope (jTsuneta et al 



upt 

3] 



2008) onboard the Hinode satellite (jKosugi et al.ll2007l ) . iLitesI (|2009l ) put forward this mechanism to 
explain how mass is discharged from magnetic field lines in emerging flux regions. In the simulation, 
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mass is subducted in a similar fashion though in the 3D case plasmoids take the place of O-loops. 



The amount of vertical unsigned flux at the surface may erode depending on the height at which 
reconnection occurs. Specifically, if reconnection happens exactly at the surface, the unsigned flux 
will immediately decrease. If reconnection occurs above r = 1, the unsigned flux will decrease 
on ly when resulting O-loops (plasmoid s) submerge. This second scenario has indeed been reported 
by llida. Yokovama. &: Ichimotol (|2010l ). If reconnection occurs below the surface, the surface un- 
signed flux will decrease if the remnant U-loop (counterpart to the O-loop after reconnection) rises 
above the surface. This is less likely since the U-loop will still be anchored to downflows. We find 
in our simulation that reconnection occurs both in the photosphere and in the convection zone. 
However, the unsigned flux (at any one horizontal level) erodes due to retraction of inverse O-loops 
(plasm oids) . The pro cess described has b een n umerically modeled by llsobe. Tripathi. fc Archontis 
(|2007l . in 2D) and by I Archontis fc Hoodl (|2009i . in 3D) for explaining the origin of Ellerman bomb s 
resulting from magnetic reconnection occurring in emerging flux regions (jPariat et al.ll2004l . |2009| ). 



The right panel of Fig. [8] shows an example of a collection of sinking U- loops in the simulation. 
The downward transport of such U-loops, which allows the discharge of m ass from the rising 
magnetic field, is rel ated to the phenomenon known as turbulent pumping (jTobias et al.l 1 19981 ; 
Brummell et al.1 l2002l ). Turbulent pumping manifests itself in terms of vertical component the 
Poynting flux of magnetic energy, which for ideal MHD is 



S z 



1 

4-7T 



[B 2 X + By) - B z (B x v x + B yVy )dxdy, 



(9) 



where the first term in the integrand represents the bodily ascent (or descent) of horizontal fields 
and the second term represents the Poynting flux due to horizontal flows. Figure [9] shows plots of 
the Poynting flux through two horizontal planes (z = and z = —4 Mm) as functions of time. The 
individual contributions from vertical and horizontal flows are also plotted. The traversal of the 
rising magnetic structure through the z = —4 Mm plane appears as a positive hump in both the 
emergence term and the combined Poynting flux between t = 1.5 and t = 4 hr. Thereafter the effect 
of turbulent pumping takes over and the net Poynting flux due to vertical flows is negative. At the 
surface (z = Mm), the convective downflows are much stronger (in terms of Mach number, they 
approach 0.1) and the effect of turbulent pumping on the Poynting flux is pronounced throughout 
the emergence and post-emergence phase. It should be noted that although the Poynting flux (a 
flux of energy) has a negative sign, it does not necessarily imply that net magnetic flux itself is 
migrating downwards. 



3.4. Growth of Magnetic Flux Content in the Developing Spots 

As already discussed in the previous section, mean horizontal outflows persist through the 
simulated duration of AR formation. Given such mean flows, how does the magnetic flux at the 
photosphere migrate inwards to yield coherent spots? To examine this, the magnetic and velocity 
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fields (for a fixed time) were decomposed into the sum of azimuthal averages and corresponding 
fluctuating components, namely B(r, 6, z) = B(r, z) + B'(r, 6, z) etc, where the bar and prime 
denote an azimuthal average and the fluctuation about the average, respectively. For ideal MHD, 
the mean field induction equation is 

— = V x (v x B) + V x £, (10) 
at 



where £ = v' x B' is the (azimuthally-averaged) mean-field electro motive force resulting from cor- 



relat ions between the magnetic field and velocity field fluctuations (jKrause fc Radlerlll980l : iRadler 



1980). To examine how (vertical) magnetic flux accumulates in the growing spots, consider a cir- 
cular area C of fixed radius R, at height z = and centered at r = 0. The time rate of change of 
the total magnetic flux crossing this circular area is 

# c - J e f dS (11) 



V x (v x B + £j dS, (12) 
vxB + fVdi, (13) 



where 



dC 

$m + $f, (14) 



$ m = 27rR[vxB] e , (15) 
= 2TrR(v z B r -B z v r ), (16) 
$ f = 2irR£ e , (17) 



= 2TrRv' z B' r - B' z v' r . (18) 

Equation (I14D indicates that the magnetic flux enclosed within C changes as a result of (1) advection 
of the mean magnetic field B by the mean flow v (first term on the r.h.s.), and (2) correlations 
between their fluctuating components. 

Figure [10] shows plots of the mean magnetic field and of the flux transport rates defined in 
Eq. (I14p . To smooth out fluctuations in time, the azimuthal means and azimuthally-fluctuating 
quantities were also averaged between t = 9.3 and t = 10 hr in time. This remains consistent 
with the derivation above since the operations of averaging temporally and azimuthally commute. 
At this stage of the simulation, the mean magnetic field within a radius of 8 Mm has already 
reached kG values and we focus our attention on how the weaker field at larger radii is transported. 
Inspection of the mean velocity field reveals that the contribution from the radially directed outflow 
(v r > 0) dominates such that magnetic flux is transported away from the spot. However, it is found 
that correlations between the velocity and magnetic fields (i.e. £ ) provide a means for the inward 
migration of vertical field. The amplitude of <3?f is in fact larger than that of 3> m and the net result 
is an accumulation of flux in the spot. 



- 11 - 



The streaming of opposi t e polarity (vertical) field in opposite horizont al directions (jStrous et al 



1996; Strous k, Zwaan 



19991 ; iBernasconi et al.l 12002 ; ICheung et al.l 120081 ) provides the systematic 



correlation needed for accumulation of magnetic flux in the spot. This begs the question of what 
is the underlying driver of the coun t er-streaming motion. The left panel of Figure [8] shows how 
the granular expulsion (jWeissI 119661 ; iHurlburt &; Toomrd Il988l ) of magnetic flux from serpentine 
field lines leads to a migration of positive polarity in one horizontal direction and negative polar- 
ity field in the opposite direction (i.e. a non-zero £). The presence of sea serpents, however, is 
an insufficient condition for flux accumulation in spots. For instance, consider a scenario wherein 
uniform magnetic field (with a net horizontal component) threads a series of granules. Although £ 
is non-zero, the translational symmetry of the setup causes the loop integral §q C £ ■ dl to vanish so 
that there is no net change of flux. 

The missing ingredient is some underlying large-scale structure which remains coherent over 
time-scales much greater than the granulation turnover times at the surface. Here it is provided 
by the coherent subsurface roots of the nascent active region (which is illustrated by the anchored 
field line in the left panel of Fig. [8]). The upper four panels of Fig. [TT] show radial profiles of 
the azimuthally-, temporally- and spatially- (over a height of 540 km about z = 0) averaged forces 
about the developing spot. For both positive and negative polarities, the gas and magnetic pressure 
gradient forces roughly (but not exactly) balance each other. The profile for the magnetic tension 
force (B • VB/47r) shows a tendency to accelerate positive (negative) polarities inward (outward). 
This trend is even clearer when one looks at the net force (-V(p gas + B 2 /8tt) + B • VB/47r) and 
is consistent with the counter-streaming behavior of the opposite polarities (see bottom panel of 
Fig. [TT]) . From this analysis, we conclude that the organization of the surface flux fragments into 
the two spots results from the presence of the coherent surface roots which influence the surface 
dynamics via the Lorentz force (especially the tension force). 



This result closely resembles the tethered-balloon model of lSpruitl (Il98ll ). In this model, emerg- 
ing magnetic loops consist of buoyant gas in the crests and neutrally-buoyant or anti-buoyant gas at 
footpoints anchored below the photosphere. Like tethered balloons, the equilibrium state is reached 
when buoyant elements hover vertically above their footpoints. The evolution of a field-line segment 
from the oblique to the vertical is mediated by the tension force. In our case, the interaction of 
granular convective flows with emerging magnetic fields, resulting in the und ulation of field lines , 
flux expulsion to inte rgranular lanes and subsequent convective intensification (jCheung et al.ll2008l : 
Danilovic et al.ll2010l ) add a layer of complexity to the problem. Nevertheless, the present analysis 
demonstrates that Spruit's model is qualitatively compatible with our simulation results. 



3.5. Light Bridge Formation and Disappearance 

As coalescence of flux concentrations progresses, ambient material with relatively weak field 
can become entrained in between adjacent high field-strength concentrations. An extended lane of 
entrained material trapped between neighboring regions of strong field appears as a light bridge 
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in intensity maps. Figure [12] shows an example from the simulation. In this case, the horizontal 
convergence of strong field concentrations (manifested as dark umbrae in continuum intensity) 
traps ambient material with weak field to form a light bridge. Three vertical cross-sectional cuts 
oriented perpendicular to the light bridge show the vertical velocity distribution and magnetic field 
strength along parts of the bridge. At its end, located deep inside the spot ([x,y] = [—18,0] Mm), 
the light bridge has a width of a few hundred km. The corresponding cross-sectional plots of v z 
and \B\ show that it is associated with an upflow of plasma with relatively weak magnetic field (a 
few hundred G at maximum). The upwelling locally lifts the T500 = 1 surface by about 300 km 
relative to the adjacent umbra. The expansion with heigh t of the neighbor ing magnetic field pinches 
the upflow to force a cusp-like structure n ear rsno = 1 (INordlundl 12004 ) . This is consistent with 
the observational study of light bridges by iLites et al.l (|2004l ). Just below T500 = 1, the pinching 
accelerates the upflow from 0.5 — 1 km s _1 to a few km s _1 . Above T500 = 1, the upflow speed 
decreases. This stagnation-like point partially traps low entropy material (radiatively cooled) and 
leads to a central dark lane withi n the light bridg e. The slower speeds at the surface are consistent 
with the observational results by iRimmeld (|1997l ). Cool material escape by feeding the downflows 
flanki ng the upflow. All these processes are similar to those occurring in simulations of umbral 



dots (ISchussler fc Vogleill200fl ) 



Vertical cross-sections across wider segments of the light bridge (Fig. fT2j) tell a similar story. 
The vertical cross-section given in the middle panel with weak field, also flanked by a pair of 
periphery downflows. At an even wider section of the bridge with a width of a few Mm, the 
convective flow consists of a few overturning cells instead of a single one. 

As the strong field concentrations continue to converge, the light bridge eventually disappears 
(at around t = 17 hrs) by means of outflows channeled toward the edge of the spot. However, 
until the time when we halted the simulation (t = 37 hrs), transient light bridges and umbral dots 
continued to appear in the pair of spots due to the tendency of entrained material to attempt to 
escape. Towards the later stages of the simulation, however, the frequency (and size) of light bridges 
diminished together with the amount of entrained material available for light bridge formation. 



Togethe r with previous models of umbral convection (jSchiissler fc Voglerl 120061 ) and sunspot 



simulations (jHeinemann et al.l 120071 ; iRempel et al.ll2009al lbl) suggest that umbral dots, penumbral 



filaments and light bridges are all manifestations of overturning magnetoconvection. 



4. Discussion 

In this paper, we presented the first radiative MHD simulation of the birth of an AR on the solar 
surface. To mimic the emergence of magnetic flux from the convection zone to the photosphere, 
a magnetic semitorus was kinematically advected upward through the bottom boundary, which is 
located at 7.5 Mm below the photospheric base. Here is a list of interesting findings regarding the 
physics of how AR formation occurs in the simulation. 
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The magnetic field strength B of the rising plasma in the emerging flux region roughly follow 
the scaling relation B oc g 1 ^ 2 , where g is the mass density. This scaling relation is consistent with 
the scenario that rising plasma preferentially expand in the horizontal directions. This horizontal 
expansion disperses field over a large horizontal area. As time progresses, the dispersed magnetic 
fragments coalesce into gradually larger concentrations (in terms of flux content). When compact 
magnetic concentrations attain a flux of ~ 10 19 Mx or more, they appear as dark pores in intensity 
images. 

In order for dispersed field to become compact again, excess mass must be removed from the 
original emerging structure. This is facilitated by convective downflows, which act on horizontal 
field to form U-loops anchored below the surface. Magnetic reconnection within these loops trans- 
po rts mass aw ay from rising field lines. This physical mechanism is the same as the one suggested 
by iLites! (120091 ) for explaining how the mass is discharged from the rising magnetic field in observed 
emerging flux regions. 

Following the removal of mass from the original magnetic structure, magnetic flux must 
migrate into two concentrated trunks to form spots. It is found that the azimuthally- averaged 
flows around the developing magnetic concentrations (corresponding to spots) contribute to the 
erosion of magnetic flux. However, correlations between the fluctuations of velocity and magnetic 
fields result in net inward migration of magnetic flux with the same sign as that of the spot. Such 
correlations result from the counter-streaming motion of opposite polarity fragments. This re- 
organization of the dispersed flux fragments is due to the presence of coherent subsurface magnetic 
roots, which influence the surface dynamics via the Lorentz force. In this s ense, our result is akin 
to the tethered-balloon model of sunspot formation suggested by ISpruitl (|198ll ) . We emphasize 
that these streaming m agnetic polarities are more akin to the moving dipolar features (MDFs, 



Bernasconi et al.l |2002| ) found in emerging flux regions (jStrous et al. 



199 



Schlichenmaier et al 



20101 ) than the so-call ed Type I moving magnetic feature s (MMFs), which likely also have serpentine 
field line geometries (jSainz Dalda &: B ellot Rubio 20 081; iKitiashvili et al.ll2010l ) but which may be 
associated with the decay of sunspots (IKubo. Shimizu. Tsunetall2007l ). 



As the dispersed magnetic field collect together to form coherent spots, plasma with relatively 
weak field (at most a few hundred G) and high entropy may be entrained between regions of 
predominantly vertical field of a few kG strength. These strong field regions correspond to dark 
umbrae in intensity maps whereas the entrained weak field plasma in between manifest itself as a 
light bridge. Cross-sectional cuts of the vertical velocity and magnetic field strength distribution 
perpendicular to a light bridge in the simulation shares strong resemblance with similar cross- 
sectional cuts through umbral dots and penumbral filaments. This suggests a common convective 
origin for all three intensity features. 

In this paper, we have focused our attention on analyzing how physical processes occurring 
in our radiative MHD simulation allow for the formation of an AR with certain observational 
properties (e.g. elongated granules and mixed polarity patterns in the emerging flux region, pore 
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formation, transient light bridges etc). In a follow-on paper, we intend to carry out more detailed 
comparisons with observations of emerging flux regions. Such an exercise would be especially 
timely considering the rise of solar activity and the availability of vector magnetograms from both 
the Helioseismic Magnetic Imager onboard the Solar Dynamics Observatory and from the Hinodc 
Solar Optical Telescope. 

This work was made possible by NASA's High-End Computing Program. The simulation 
presented in this paper was carried out on the Pleiades cluster at the Ames Research Center. 
We thank the Advanced Supercomputing Division staff for their technical support. This work was 
supported by NASA contract NNM07AA01C at LMSAL. M. Rempel is partially supported through 
NASA grant NNH09AK02I to the National Center for Atmospheric Research. NCAR is sponsored 
by the National Science Foundation. 
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Fig. 1. — Time sequence of continuum intensity images at 500 nm (left) and synthetic longitudinal 
magnetograms (sampled at tr ss = 0.1, right). The images show the full horizontal extent (92 x 49 
Mm) of the simulation domain. 
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Fig. 2. — The upper panel shows a synthetic magnetogram (sampled at tr = 0.1) of the simulated 
emerging flux region at t = 6.6 hrs. The lower panel shows the vertical cross-section of the magnetic 
field strength (log 10 B) along y = (delineated by red dashed line). 




Fig. 3. — Same as Fig [2] but at t = 15.3 hrs. The two polarities of the active region at this time 
are much more compact. 



(a) 



(b) 



Fig. 4. — Magnetic structure of the (semi)torus as it erupts onto the solar surface (t = 4 hrs). 
The same set of field lines are shown in both panels. The left panel shows the distribution of the 
vertical component of the magnetic field (scaled between ±17 kG) at z = —7.5 Mm. At this depth 
the opposite polarities are in the form of two coherent roots. The panel on the right shows the 
corresponding magnetic map at z = (scaled between ±2 kG), which highlights the serpentine 
nature of the magnetic field lines near the surface as a consequence of the interaction between the 
emerging field and the granular convective flows. 




Fig. 5. — Velocity field below the surface (z = — 2 Mm) as the emerging magnetic field is traversing 
this layer (t = 2.8 hrs). Red and blue colors indicate down- and upflows, respectively, and the 
arrows indicate the horizontal velocity. The outflows at the periphery of the emerging flux region 
reaches speeds of up to ~ 4 km s _1 whereas upflow velocities are closer to ~ 0.5 km s _1 . 
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Fig. 6. — Joint probability density function of the horizontal field strength. The values were 
sampled at the x = vertical plane between t = 1 and t = 8 hours. The solid and dashed lines 
respectively show the scaling relations B oc g 1 / 2 and B oc g 2 / 3 with Bo = 21 kG and go = 4.2 x 10~ 4 
g cm -3 . 
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Fig. 7. — The emerging flux tube leads to a enhancement of pressure which drives outflow away 
from the emergence site. The three panels in this figure show, respectively, the time evolution of 
the surface (photospheric base, z = 0) magnetic field strength, relative gas pressure perturbation, 
and x— component of velocity averaged over the band y G [—2, 2] Mm. In all three panels, the green 
contours indicate enhancements of the gas pressure (relative to (-P ga s) = 9 x 10 dyne cm -2 ) by 
25% and 50% (for the purpose of having fewer contours, the pressure values have been smoothed 
in time with a Gaussian filter with a = 30 min). 
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(a) (b) 

Fig. 8. — Mechanism for removal of mass and unsigned flux from the surface, (a) Schematic 
representation of how mass is removed from emerging magnetic field lines in a 2D scenario. In 
addition to undulating field lines, convective flows expel emerged flux (indicated by ovals labeled 
with positive and negative signs) from granular upflows. (b) 3D rendering of near-surface field lines 
in the simulated emerging flux region. Field lines in the upper panel are colored in accordance 
with the local density perturbation (about horizontal mean) with dark blue indicating density 
enhancement. Field lines in the lower panel are colored according to the vertical component of the 
momentum with red indicating downflowing material. 
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Fig. 9. — Poynting fluxes of magnetic energy through horizontal planes at z = and z = —4 
Mm. In both plots, the dashed curves and diamonds indicate the contributions by the first (bodily 
emergence or submergence of horizontal fields) and second (shearing by horizontal flows) terms in 
Eq. (|9|), respectively. The solid curves indicate the sum of the two contributions. 
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Fig. 10. — The top panel shows the azimuthally- and temporally-averaged (between t = 9.3 and 
t = 10 hr, sampled at z = 0) magnetic field strengths (solid and dashed lines indicate the vertical 
and radial components, respectively) as functions of the radial distance from the axis of the positive- 
polarity spot in the simulation. The lower panels shows the corresponding quantities as defined in 
Eq. (fT4|) : <& m (dotted) , $f (dashed) and their sum (solid) . Positive values of <3? indicate an increase 
of flux (within the area r < R) with time. 
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Fig. 11. — The panels above show the azimuthally-, temporally- and spatially- (over a height-range 
of 540 km about z = 0) averaged profiles of the radial components of the magnetic pressure gradient 
force (— VB 2 /8ir), gas pressure gradient force (— Vp gas ), magnetic tension (B • VB/47r), net force 
and velocity (v) for gas with positive polarity field (B z > 0, thick solid lines) and negative polarity 
field (B z < 0, thick dashed lines). The developing spot centered at r = has positive polarity. 
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Fig. 12. — Structure of a light bridge - The left panel shows the continuum brightness of one of 
the spots in the simulated active region (t = 15.4 hrs). The formation of this light bridge is due 
to buoyant material entrained between adjacent regions of high magnetic field strength. Vertical 
cross-sections of the vertical velocity and magnetic field strength are shown for three cuts across 
segments of the light bridge of varying width. The green lines in the cross-sectional panels indicate 
the T500 = 1 surface. 



